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Abstract 

We analytically derive the spectrum of gravitational waves due to magneto- hydrody- 
(SJ namical turbulence generated by bubble collisions in a first-order phase transition. In 

contrast to previous studies, we take into account the fact that turbulence and magnetic 
fields act as sources of gravitational waves for many Hubble times after the phase tran- 
sition is completed. This modifies the gravitational wave spectrum at large scales. We 
also model the initial stirring phase preceding the Kolmogorov cascade, while earlier 
works assume that the Kolmogorov spectrum sets in instantaneously. The continuity 
in time of the source is relevant for a correct determination of the peak position of 
. ^ the gravitational wave spectrum. We discuss how the results depend on assumptions 

^ about the unequal-time correlation of the source and motivate a realistic choice for 

it. Our treatment gives a similar peak frequency as previous analyses but the ampli- 
tude of the signal is reduced due to the use of a more realistic power spectrum for 
the magneto-hydrodynamical turbulence. For a strongly first-order electroweak phase 
transition, the signal is observable with the space interferometer LISA. 
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1 Introduction 

Cosmological observations are often a search for 'relics' of the early universe. Relics allow 
us to infer the physics at a time when the Universe was much hotter and much denser than 
today. A famous example of this is the cosmic microwave background (CMB) [1], which 
literally represents a photograph of the Universe at the time when CMB photons decoupled. 
Another very promising relic which has not yet been observed is a gravitational wave (GW) 
background. Since GWs interact so little with matter and radiation, they propagate freely 
immediately after generation and therefore allow us a direct observation of the Universe at 
the time of their production. GW backgrounds have been proposed from inflation [2], from 
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braneworlds [3] , from topological defects [4] , from reheating after inflation [5] and from first 
order phase transitions [6,7]. In this last case, which is the topic of this paper, at least three 
different sources of GWs have been identified: the collisions of broken phase bubbles [8-14] , 
fluid turbulence [15-20] and magnetic fields [21-24]. In this paper we concentrate on the 
latter two aspects which cannot be truly separated since the magnetic fields are processed 
and amplified by the turbulent fluid flow. 

The GW signal from a first-order phase transition has a characteristic frequency of the 
order u = ck^a^ja^ = e^-f/^a*/ 'a where e = L^H^/c is the size of the largest bubbles in 
units of the horizon size at the transition, 1/H*. Its value depends on the particle physics 
model but for a strong first-order phase transition, we will set e ~ 0.01. Quite remarkably, for 
a potentially first-order electroweak (EW) phase transition at T* ~ 100 GeV, this frequency 
is around a MilliHertz which is the frequency of best sensitivity of the planned GW satellite 
LISA (Laser Interferometer Space Antenna) [25], meaning that LISA is potentially a window 
on EW and TeV scale particle physics [26-28]. 

The main difference of our approach with earlier works is that we consider a long-lasting, 
magneto-hydrodynamic (MHD) turbulent source. The collision of bubbles of the broken 
phase causes an injection of energy in the primordial fluid. Since the kinetic Reynolds 
number of the fluid is huge, for instance ~ 10 13 at the EW epoch (c.f. section 3.5), turbulent 
motion sets in rapidly in the fluid. The magnetic Reynolds number being also very large, this 
leads to the amplification of magnetic fields generated during the phase transition [29] , and 
MHD turbulence develops. Once the phase transition is over the source of energy injection 
stops. However, MHD turbulence does not cease immediately but decays like a power law. In 
previous studies [15-18,20], the free decay of the turbulent velocity power spectrum has been 
ignored: it was assumed that turbulence was active only during the completion of the phase 
transition. However, since the value of the Reynolds number is very high, the dissipation is 
not sudden and the fluid can remain turbulent during many Hubble times. 

The power law decay in time of MHD turbulence is well established theoretically [30-35] , 
experimentally [36] and by numerical simulations [37-39]. However, there is no consensus 
on the actual value of the power law exponent: Kolmogorov theory predicts a faster decay 
than what observed in general in numerical simulations. On the other hand, both analytical 
analyses and numerical simulations do agree on one point: that, in the absence of helicity, 
the large scale part of the MHD turbulent spectrum is constant in time. Without inverse 
cascade the energy is dissipated on the very small scales, the correlation length grows in 
time, but wavenumbers much smaller than the inverse correlation length are not affected by 
the evolution. This is important for MHD turbulence in the early universe. The radiation 
dominated universe is characterised by a finite causal horizon, beyond which the turbulent 
motions cannot be causally connected. As previously demonstrated [17,40], the existence of 
this causal horizon in a cosmological setting implies that the real-space correlation function 
of the stochastic velocity field has compact support. The fact that the real-space correlation 
function necessarily vanishes at large scales, together with the property of divergence freeness 
satisfied by the incompressible turbulent flow and by the magnetic field, entails the formation 
of a Batchelor spectrum for the turbulent velocity field and the magnetic field [17]. Namely, 
taking for example the velocity field, the large scale part of the spectrum grows as P v (k — > 
0) oc Ik 2 , where / ~ (v 2 )L 5 is the Loitsyansky's integral, (v 2 ) being the typical velocity of 
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the largest eddies and L their size {i.e. the largest scale on which turbulence develops, the 
correlation scale corresponding to the bubble diameter in our context). Given this form of 
the large scale part of the power spectrum, if it has to be constant as predicted by the theory 
and observed in numerical simulations, then necessarily / must be constant in time (which 
is also required by the Navier-Stokes equation) [41,42]. 

In the following, we will assume that the kinetic energy (v 2 ) and the correlation length L 
evolve in such a way, as to maintain the product / ~ (v 2 )L 5 constant. We define the power 
law exponent 7 such that L ~ t 7 and (v 2 ) ~ t~ 57 . For generality, the exponent 7 is kept 
unspecified in the analytical formulae, but for the numerical results we substitute the value 
7 = —2/7. According to Kolmogorov theory, in fact, the constancy of / together with the 
energy decay equation d(v 2 )/dt ~ —(v 3 )/L lead to the Kolmogorov decay laws: the decay 
of the kinetic energy with time as (v 2 ) ~ t~ 10//7 and the growth of the correlation scale as 
L ~ t 2 ! 1 (see for example [30]). As explained in section 3.3 in the following we also assume 
equipartition between the turbulent and magnetic energy densities (v 2 ) ~ (b 2 ): consequently, 
we assume the same decay law also for the magnetic field, i.e. (b 2 )L 5 = constant. 

The value 7 = 2/7 that we use in the numerical estimates has also been derived on 
the basis of detailed theoretical arguments, as for example in [35]. On the other hand, 
as previously mentioned, numerical simulations observe a slower decay for the kinetic and 
magnetic energies, close to (v 2 ) ~ (b 2 ) oc t^ 1 [37-39]. However, it is not clear whether 
numerical simulations can efficiently model the conditions of MHD turbulence in the early 
universe, which is characterised by the presence of a causal horizon and develops at extremely 
high Reynolds number, of the order of 10 13 (c.f. the discussion at the end of section 3.4). 
Therefore, in our analysis we have chosen to follow the theoretical picture of Ref. [35], 
which in addition leads to a conservative estimate of the production of G Ws: in fact, MHD 
turbulence which decays slower would be active as a source of G Ws for a longer time. 

In the following we assume this model of free decay for the turbulence, but to this 
'absolute' time behaviour we also add the exponential de-correlation proposed in [43], to 
express the time de- correlation of the velocity field on a given scale as time goes by. We take 
the characteristic de-correlation frequency on a given scale to be the eddy turnover time at 
that scale, and assume a Gaussian functional form to express the dependence of the power 
spectrum on time difference. 

Another new point of our analysis is that the sources of GWs, turbulent kinetic energy 
and magnetic field, are continuous in time. The importance of having a continuous source 
and its consequences on the position of the peak of the GW spectrum have been analyzed 
in [14] for a short lasting source: the collision of bubbles. Here, we analyze also the long 
lasting, MHD turbulent source. In order to do so, we need to modify the Kolmogorov decay 
laws, and insert an initial phase in which the proper turbulent cascade has not yet begun: 
during this phase, we assume that the kinetic energy starts from zero and grows linearly in 
time, up to when stirring is over and the free decay of turbulence starts. The linear increase 
has been observed in MHD simulations [44] , and is also satisfied in simulations of the bubble 
collision source (c.f. [13, 14]). We assume that the evolution law of the stirring scale L 
remains equal to the Kolmogorov decay law also during the initial phase. 

Contrary to previous analyses, we also model the MHD turbulence spectrum using a for- 
mula which smoothly interpolates between the large scale behavior, determined by causality, 
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and the small scale one, given by the MHD cascade (see also [24]). This model gives a more 
realistic estimate of the amplitude of the MHD spectrum at the peak, which in turns deter- 
mines the final amplitude of the GW spectrum. With this improved MHD spectrum, the 
GW peak amplitude is more than one order of magnitude smaller than previous estimates. 

The paper is organized as follows. In the next section we discuss a toy model for the 
source to illustrate the difference between a short-lasting and a long-lasting source of GWs. In 
Section [3] we discuss the properties of the MHD turbulence, and we determine the anisotropic 
stress power spectrum from this source in Section |4j We then define the time after which we 
may neglect the GW source and determine the final GW spectrum in Section [5] We discuss 
our results in Section [6] and conclude in Section [7j A discussion of the fluid viscosity and 
some technicalities are given in appendices for completeness. 

Notation: Unless otherwise stated, we use comoving variables: t denotes conformal time, 
the energy injection scale L, the Kolmogoroff microscale A (the endpoint of the Kolmogoroff 
spectrum), and k are respectively comoving distances and wavenumber. The index * indicates 
the time of the phase transition, while the index o indicates today. We normalize the scale 
factor a(to) = 1. TC denotes the conformal Hubble parameter, Hq = /iol00km/s/Mpc is 
the Hubble parameter today, and the critical energy density is p c = p c (t ). The radiation 
energy density parameter today is /ig^rad,o — 4.2 x 10~ 5 . This value includes three types 
of neutrinos. As we shall see, this is the relevant quantity since neutrinos (with standard 
masses) are still relativistic at matter-radiation equality. 

2 A stochastic gravitational wave background of cos- 
mological origin 

We consider a Friedmann universe with flat spatial sections. The tensor metric perturbations 
are defined by 

ds 2 = a 2 (t)[-dt 2 + (Sij + 2h ij )dx i dx j ] . (1) 

In this work we want to determine the GW energy density power spectrum given by (see 
e.g. [21]) 

dfl GW k 3 \h\ 2 k 5 \h'\ 2 



dlogk 2(27r)3Gp c a 2 2{2TxfGp c a? ' Pc pM ' (2) 
where ' = and ' — 4-, x — kt, and the GW energy power spectrum is defined as 



(MM)^(q,t)> = (27r) 3 <5(k-q)KM)| 2 . (3) 

Here (■ ■ ■) is an ensemble average over the stochastic process which generates the GWs. The 
Dirac delta function 5(k — q) is a consequence of statistical homogeneity. 

Once the source has decayed and the wavelength under consideration is inside the horizon, 
the GW energy density simply scales like a~ A . Hence the GW energy spectrum scaled to 
today becomes 



dQ 



GW 



dlog k 



dQ 



GW 



dlog k 
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To evaluate the GWs emitted by turbulent motion in the primordial fluid and by a 
magnetic field we need to determine the tensor-type anisotropic stresses of these sources. 
They source the evolution equation for the GW perturbations, 



hij + 2Hh 



+ eh i:i = 8nGa 2 T^ T) (k,t) 



(5) 



In this section we consider in all generality a relativistic source, and we solve the wave 
equation in two cases: a long lasting source (i.e. many Hubble times), and a short lasting 
one (i.e. significantly less than one Hubble time). We introduce the transverse traceless 
tensor part of the energy momentum tensor of the source as 



T^ T) (k, t) = (p + p)U lj (k, t) so that 



87rGa 2 T 4 f T) (k, t) = m 2 flij(k, t) , (6) 



"ij ki kj > 



where we denote the dimensionless energy momentum tensor with a tilde: 11^ (k, t) 
(P il P jm -l/2P ij P lm )T lm (k, t). The projection tensor 

1/2^™ with P<. 

projects onto the transverse traceless part of the stress tensor, n includes any time depen- 
dence other than the basic radiation-like evolution. We assume that the source is active only 
during the radiation-dominated era, where p = p/3. During adiabatic expansion g(Ta) 3 = 
constant so that 



P(t) 



Prad.O / go 



a 4 (t) \g(t) 



1/3 



and 



(7) 



where g(t) is the number of relativistic degrees of freedom at time t. 



2.1 Long-lasting source 

Let us first concentrate on the more general case of a long lasting source. To solve Eq. ^ 
we set 7i = 1/t, neglecting changes in the number of effective relativistic degrees of freedom. 
In terms of the dimensionless variable x = kt Eq. ^ then becomes 

hf: j + 2 h f + h i j = ^U ij . (8) 

We consider a source that is active from time t m to time £g n , which in the long lasting case 
can span a period of many Hubble times. For t > ifi n , we match the solution of the above 
equation to the homogeneous solution, n^- = 0. Assuming further that we are only interested 
in modes well inside the horizon today, x ^> 1, the resulting GW energy power spectrum 
becomes 

i,//, \ t 2 8 f Xfin dxi f Xfin dx 2 , ,~ /7 . . * 

\h (k, X > X fin ) / / COS{X2 — Xi)U.(k,X 1 ,X2) 2>1, (9) 

x 2 J Xin Xl J Xin x 2 

x\ = kti, X2 = kt%, and Ti(k, X\, x 2 ) denotes the unequal time correlator of the source, 

(n^MOffefota)) = (2n) 3 5(k-q)fl(k,kt 1 ,kt 2 ). (10) 
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With Eq. Q, the power spectrum of the GW energy density parameter for a long-lasting 
source which is active between t m and tg n in the radiation era is then given by 

4 qad.O [90 \ V3 ,3 r n dXi f fi " dx 2 ~ 

= ——— [ — k / / cos(x 2 - xi)U{k,x 1 ,X2) ■ (11) 

3vr 2 \g &n J J x . m xi J x . n x 2 

This result is completely general for modes well inside the horizon today; it reduces the 
computation of the GW spectrum to the determination of the unequal-time correlator of the 
tensor-type anisotropic stress, fl(k, Xi, x 2 ). 



dVt GW 
dlogk 



2.2 Short-lasting source 



If the source is active for a short interval of time, essentially only during the phase transition, 
one can neglect the expansion of the universe during the time of action of the source, and 
match the solution so obtained with the one of the homogeneous equation in which expansion 
is taken into account. For this kind of source, we set tg n = ti n + At, with At/t in <C 1. Solving 
the wave equation (|5j without expansion term, amounts to neglect the time-dependence of 
the factors \jx\ and \jx 2 in Eq. ^ or (11) during the active period, so that 



\h'(k,x > x &n )\' 



fi ii 



dx\ 



''in 



fin 



dx 2 cos(x 2 — xi)il(k, xi,x 2 ) 



(12) 



Also this solution applies for modes inside the horizon today. The energy spectrum now 
becomes 



dVL 



GW 



d\og k 



4 Q 



rad,0 



3vr 2 



9o 

9&n 



1/3 



nlk 



■'fin 



dx\ 



■' tin 



dx 2 cos(x 2 — Xi)H(k, Xi, x 2 ) ■ (13) 



Obviously, the general long- lasting case reduces to this result if (tfi n — t in )/t in <C 1. Summarizing: 



long-lasting source {e.g. MHD turbulence) 

(14) 

(?) 3 k J2 n dx ± JZ n dx * co < x * ~ Xl ) Xl ' X2 ) 

short-lasting source (e.g. bubble collisions) 
with x\ = kti, x 2 = kt 2 , x in = kt in , Xfi n = kt^ n and A = g^^rad.o^o • 

2.3 Solutions for a simple source 

In this section we analyze the difference between the GW spectrum generated by a short 
lasting and a long lasting source in a simple example which can be treated analytically. We 



dQcwh-l 



dlog k 



Al 
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find analytical solutions for the GW energy density spectrum (14). The general behaviour 
of the solutions in this simple case is illuminating, as it is similar to the case of the evolving 
source that we will treat in the rest of the paper (MHD turbulence). 

We consider a tensor source IT with an equal time power spectrum which depends on time 
solely via a function modeling the turning on and off of the source: H(k,t,t) = U(k)f 2 (t). 
For the equal time power spectrum of the tensor source, we take a form which is motivated 
by turbulence and magnetic fields, see Sec. |4j 



rad 



2 



U(K,t,t)=[^ \ L 3 S(K)f(t), (15) 



where Q$ denotes the (radiation like) energy density of the source normalised to the critical 
energy density today (so that the ratio f2,s/f2 rac i is time- independent), K = Lk/2TT is a 
dimensionless wavenumber, L is a characteristic scale of the problem, and S (K) models the 
scale dependence of the source. The continuous function f(t) vanishes at both, t in and t nn 
and describes the switching on and off of the source. The source is active during the time 
interval t fin — t in , which can be long or short compared to the initial Hubble time, TC^ = t in . 

As shown in Ref. [14], the time continuity in switching the source on and off can be 
relevant for the GW spectrum. Inspired by results from numerical simulations of bubble 
nucleation during a first order phase transition [10,13], we choose the function f(t) to be 
continuous but not differentiable at the initial and final times. The effect of this choice on 
GW spectra for short duration sources is discussed in Ref. [14]. Here we shall also study its 
effect on sources of long duration. We set (see Fig. [T| 



/(*) 



if t < t 

2(t-t in ) 
At 



if t in < t < t in + At/2 
1 if ^4- At/2 < t < t^- At/2 (16) 

if t^- At/2 
if t > t 



fin 



For a short lasting source, At is also the duration of the source (c.f. Sec. 2.2 and Ref. [12]), 
while a long lasting source is typically active for a much longer period of time, and At is the 
characteristic time of turning on and off. As we show below, the effect of introducing time 
continuity on the resulting GW spectrum is relevant only in the 'coherent' case (see Fig. [3J. 
In this case, the power at small scales is less than for a discontinuous source. 



To compute the solution according to Eq. (14), we need the source power spectrum at 
unequal times n(fc,ti,t 2 ). In this section, we consider two different possibilities for the 
unequal time correlators, the incoherent and the totally coherent approximations, which we 
discussed also in the case of GWs generated by bubble collisions [12]. In Sec. [3] we will 
compare these approximations with the 'top hat' ansatz [12], which turns out to be more 
realistic for the case of MHD turbulence. 

• Incoherent approximation: the source is correlated only for t\ ~ t 2 

U(K, t u t 2 ) = U(K, tutJSfa - t 2 ) At f 2 (h) . (17) 
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Figure 1: The function modeling the time dependence of the source used in Section 2.3, Eq. (16). 



Here At plays also the role of a very short characteristic time, over which the source remains 
'coherent'. With this ansatz the time integration in (14) is simple, and we obtain the 
following result for the GW energy spectrum: 



where 



dlogk 



F(t in , t nn , At) 



3tt 2 



(— ) 



n 



rati 



3 At 

tin 



K 3 S(K)F(t in ,U n ,At) 



(9g\f (2?r) 2 (At\ 2 
KgJ 3 \t in J 



long-lasting, 



short-lasting. 



(19) 



The result for the long-lasting case is expanded using At <C ti n , t nn . The full expression is 
given in Appendix |Al where we collect all analytical expressions for the convenience of the 



reader, see Eq. (97) 



In the incoherent approximation, the function F(t- m , tg n , At) resulting from the convo- 



lution of the Green function with the source, Eq. (14), does not depend on wave- number. 



The GW power spectrum is therefore simply the one of the source S(K), multiplied by the 
phase-space volume K 3 , both for long lasting and for short lasting sources. The peak of the 
GW spectrum then coincides with the one of the source spectrum. The time integration 
results in the ratio between the brief coherence time At and the initial horizon time t in . In 
the long lasting case, this factor might be close to one while in the short lasting case it is 
always much smaller. In addition, the GW amplitude of the short lasting case is suppressed 
by one more factor At /tin with respect to the long lasting one. 

• Coherent approximation: the source is perfectly correlated at all times t% and £2 

n(tf,t ls t 2 ) = vn(^0i(i^). (20) 



In this case as well, the time integration in Eq. (14) can be performed explicitely, and we 



obtain the GW energy density power spectrum (we remind that X[ n = kt- m , Xf\ n = kt nn , and 

Ax Xf[ n X[ n ) 



dlog k 



4^rad,0^0 ^ ^' 



37T 2 



n 



ra< 1 



K 3 S(K)F(x in , x fin , Ax) 



(21) 



where 



Si(x in )) 2 ] long-lasting 

(22) 
short-lasting. 

In the long lasting case we have again expanded to lowest order in Ax/x in and Ax/x Rn . Ci 
and Si denote the integral cosine and sine functions [45]. The full expression is given in 
Appendix [A} 

Contrary to the incoherent case, here the function F(x- m , xg n , Ax) depends on wave- 
number, and the resulting GW spectrum is therefore modified with respect to the one of 
the source. We plot F(x- m , x^ n , Ax) in Fig. |2j In both the long and short lasting case, 
F(xi n , Xfi n , Ax) tends to a /c-independent value for wave-numbers such that x- in <C 1. The 
constant is log 2 (xfi n /xi n ) in the long lasting case, and 7r 2 (Ax/xi n ) 2 in the short lasting one. 
In the short lasting case F remains constant up to Ax ~ 1 where it starts oscillating 
and decaying like i/(x in Ax) 2 . The long lasting case instead depends also on Xfi n : when 
Xfi n becomes larger than one, the slope changes from the constant to a mild logarithmic 
dependence on k as log 2 Xi n . Then for x in > 1 we have a decay like l/x 2 n , up to Ax > 1 
where F also starts oscillating and decaying, with a smaller amplitude than the short lasting 
case, see Fig. [2} Hence, wavelengths which are larger than the typical switching on time are 
amplified by an additional factor min{(Ax) -2 , (x; n /Ax) 2 } in the long lasting case. On the 
other hand, wavelengths smaller than the the typical switching on time are suppressed in 
the long lasting case due to the presence of interferences suppressing the signal. 

The functions F(x- m , xa n , Ax) for both the short and long lasting incoherent and coherent 
cases are shown in Fig. [2} The GW amplitude in the short lasting case is suppressed for low 
and intermediate values of k with respect to the long lasting one, both in the incoherent and 
in the coherent cases. However, while in the incoherent case this suppression is maintained 
for all k, in the coherent one interferences suppress the amplitude of the long lasting case 
for frequencies larger than the typical switching on frequency Ax > 1. 

Fig. [3] shows the effect of introducing continuity in the process of turning on and off 
a long lasting source. The conclusions drawn in the analysis of Ref. [14] are not modified 
by the long duration of the source. The incoherent spectrum is the same for a continuous 
and a discontinuous source, while in the coherent case there is a difference for frequencies 
higher than k > 1/ At corresponding to the characteristic time-scale of the source. In the 
continuous case the slope changes, becoming steeper by a factor k~ x if the turning on process 
has a kink: jit) is continuous but not differentiable. In the discontinuous case, on the other 
hand, the change of slope is absent. 

Summarizing, two general features can be deduced from this analysis. First, the GW 
energy spectrum at large scales is proportional to the phase space volume K 3 times the 
source power spectrum S(K). In the incoherent case no other wavenumber-dependence 
intervenes; in the coherent case, instead, the wavenumber at which this behaviour changes 
depends on the duration of the source, whether it is short or long lasting. Second, the 
GW amplitude in the short lasting case is smaller than the one in the long lasting case. In 
the incoherent case, this is always true; in the coherent case, this suppression can be very 
significant on large scales, however, it turns into an amplification for scales smaller than At. 



F(^Xi ny Xfi n , Ax) 



(ja_)3 [( C i(x fin ) - Ci(x in )) 2 + (Si(x fin ) - 

/go\§ 64(2vr) 2 sin 4 ^^-^)/^ 

Vg*/ X? (Zfin-Zin) 2 
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Figure 2: Comparison between a short and long-lasting source. We plot the function 
F(xi n , Xfi n , Ax) defined in Eq. (19) (incoherent case) and Eq. (22) (coherent case) as a function 
of 

%m — kt[ n . Blue, solid: long lasting coherent and incoherent cases with tfi n /ti n — 100 and 
At = 0.01 tj n . Red, dashed: short lasting coherent and incoherent cases with t nn /ti n = 1.01, 
At = 0.01 t; n . The horizontal lines correspond to the incoherent case. 




Figure 3: For a long-lasting source, comparison between the continuous (solid) and discontinuous 
(dashed) cases. We plot the function F(x- m , x nn , Ax) as a function of x ia = kt- m . The incoherent 
case is not affected by continuity as shown by the red horizontal line (the solid and dashed lines 
are superimposed). In the coherent case shown in blue, the slope changes in the continuous case 
for frequencies k > At -1 . The values of the parameters are the same as in Fig. M 



On super-horizon scales, the difference in the amplitude between the short lasting and long 
lasting sources is typically of the order of At/t in in the incoherent case and (At/t in ) 2 in the 
coherent one. 
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Examples of well-motivated long-lasting sources are turbulent fluid flows initiated by 
instabilities generated by bubble collisions. In the following, we concentrate on this particular 
source. We also study the case of magnetic fields, which represent another long-lasting source 
of GW typically expected from first order phase transitions [29]. 

3 Turbulent magneto- hydro dynamics as a source of 
gravitational waves 

3.1 General considerations 

Turbulence develops if a fluid with sufficiently high Reynolds number is perturbed. If the 
fluid is stirred on a characteristic scale L p (the subscript indicates that we use the physical 
length, not comoving length here), the Reynolds number of the flow is defined by 

Re (L p ) = ^ (23) 

where vi is the characteristic velocity on the energy injection scale L p , and v is the kinetic 
viscosity of the fluid. If Re (L p ) ^> 1, the stirring develops turbulent motions. In a first 
order phase transition, the source of stirring is bubble collision. Therefore, the characteristic 
scale of the stirring is given initially by the typical bubble size towards the end of the 
phase transition, ~ 2t> b( 5 _1 where Vb is the bubble wall velocity and the (comoving) 
duration of the phase transition (see for instance Section 4 of [12] for a detailed definition). 
The initial energy injection scale is the scale at which the largest turbulent eddies develop, 
and corresponds to the peak of the turbulent velocity power spectrum. 

The dynamics of bubble growth at late times, towards the end of the phase transition, 
allows us to determine the order of magnitude of the kinetic energy involved in the turbulent 
flow. The bubble wall can be treated as a discontinuity, i. e. a combustion front across which 
energy and momentum are conserved [46]. At the front, the velocity of the fluid in the 
rest frame of the bubble center is given by Vf = {vi — t>2)/(l — ^1^2)? where Vi and v 2 are 
respectively the incoming and outgoing speed of the fluid in the rest frame of the front (see 
e.g. Ref. [12] for more details). We assume that the typical value of the turbulent fluid 
velocity is given by Vf. Therefore, the kinetic energy of the turbulent flow is 

(v 2 ) 

Pkin = (p + p) ±y with (v 2 ) ~ v) . (24) 

Even though the velocities involved are large in the case of interest, we are using the formal- 
ism of non-relativistic MHD turbulence. Since the corresponding values of 7 are typically 
of order one, we expect that the error introduced is within the uncertainty of our calcula- 
tion. Furthermore, while the validity of the theory of non-relativistic turbulence may be 
questionable if high speeds are involved [47], it was shown in [48-50] that the Kolmogorov 
spectrum is recovered even in the relativistic cas Still, we impose for the fluid velocity 

1 It is remarkable that recent simulations of quark gluon plasma instabilities in the process of thcrmaliza- 
tion in heavy ion collisions show similarities with a Kolmogorov scaling [51-53]. 
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(v 2 ) < c 2 , where c s = 1/a/3 is the sound speed in the relativistic fluid. The fluid velocity Vf 
is completely specified once Vb is known and the ratio a = p vac / Prad,* is fixed: solving for the 
hydrodynamical equation allows one to relate the fluid velocity to the bubble wall velocity 
(see for instance [54] for more details). If the phase transition proceeds as a detonation 
(deflagration), then v\ = Vb (t> 2 = ^&)- For this paper, we choose the fluid velocity Vf = c s 
corresponding to a = 1/3 and either Vb — 0.87 (detonation) or Vb — c s (deflagration) (for 
a > 1/3 there is no deflagration solution). When we need to specify Vb, in the figures and the 
numerical values, we always consider the detonation case Vb = 0.87. It is straight forward to 
re-scale the results to lower bubble and fluid velocities. 

If the fluid is stirred on the scale L*, turbulent motions develop within a time interval 
of the order of the eddy turnover time T£. This is the characteristic time for the cascade to 
set in. Given the typical value of the turbulent fluid velocity v /, the eddy turnover time on 
the stirring scale L* is defined simply as r L ~ L*/{2vf). Since the fluid velocity is always 
smaller than the bubble wall velocity Vf < Vb, the eddy turnover time is always larger than 
the duration of the phase transition: tl > (3~ l . In the following we identify the time interval 



At given in Sec. |2.3| as At = tl- For short-lasting turbulence, this means that the source 
lasts for only one eddy turnover time. In the long-lasting case (which is the relevant one as 
we will see), this means that turbulence is 'turned on' in one eddy turnover time r^. 

In the cosmological context the fluid is ionized and has not only a very high kinetic 



Reynolds number (which we evaluate in Section 3.5 ) but also a very high magnetic Reynolds 
number R m , defined by (see Appendix |b|) 



R m (L p ) = p L , where p = - — (25) 
/i Ana 

is the magnetic diffusivity and a denotes the conductivity. High values of the magnetic 
Reynolds number require an MHD treatment of the cosmic plasma. Moreover, the magnetic 
Prandl number 

„ K m (L) v , . 

p » = iRzK < 26 > 

is much larger than unity, as we calculate in Appendix [Bj Therefore, the characteristics of 
the plasma in the early universe entail the formation of MHD turbulence. The seed magnetic 
field can be generated by several mechanisms [29], and is then amplified by the currents due 
to the turbulent flow of charged particles [31,55]. The magnetic field itself is also a source of 
GWs. In the following we make the simplifying but reasonable assumption of equipartition: 
the total kinetic energy in the turbulent motion is equal to the magnetic field energy. Note 
however, that this assumption need not hold for each wave number k, so that we can allow 
for different spectra for the magnetic field energy and the turbulent kinetic energy at small 
scales. 

In the remainder of this section we present our model for the power spectra of turbulence 
and magnetic fields, as well as their time evolution. We evaluate the Reynold number in the 
early universe and describe how it evolves with time, which will help us to determine when 
MHD turbulence is expected to end. 
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3.2 The turbulent velocity power spectrum 

The power spectrum of the turbulent velocity field at equal times is of the forirj^] 

(v t (k,t)v*(q,t)) = (2ir) 3 6(k-q)P ij P v (k,t), (27) 
where the projector Pj,- = 8^ — kikj comes from the fact that the turbulent velocity field 



v(x, t) is divergence free. An ansatz for the unequal time correlator will be given in Sec. 3.6 
In previous works, the turbulent velocity power spectrum P v (k,t) was either assumed to be 
given only by the inertial range k~ u ^ 3 [15,16,18,20], or naively determined by intersecting 
the k 2 behaviour at very small scales with the inertial range fc~ 11//3 behaviour [17]. Both 
approaches overestimate the peak amplitude of the turbulent source and thus overestimate 
the GW amplitude. In the present study, we use a more realistic, smooth function to describe 
the spectrum, proposed by Von Karman [56] (see also page 244 of [57]). It is given by the 
following interpolating formula (here in terms of comoving quantities) : 

W = C^) I? x { I £ ^ * (28) 

where we again use the dimensionless variable 

K = kL/2n (29) 
and A denotes the Kolmogorov microscale, beyond which turbulent motions are absent and 



we set the spectrum to zero. The kinetic energy of the turbulent flow is given in Eq. (24): 
using this definition, we rewrite (v 2 ) in terms of the ratio of the total kinetic energy to the 
radiation energy density, 

<* 2 > = |^. (30) 
In our numerical estimates, we will use (v 2 ) = 1/3, thus corresponding to flr/^rad = 2/9. 



In equation (28), the constant C v = 10g 5 ^ 3/2 j^yf) ~ 0.0385 comes from the normalization of 
the kinetic energy spectrum, E(k) = k 2 P v (k)/(2ir 2 ), 

«*l = £5. = ? = r dkE ( k ) = -L r dkk 2 p v {k) . (si) 

p + p 2 4Q rad J y } 2Wo 

Expression (28) smoothly interpolates between the large scale, K 2 behaviour, and the inertial 
range i^~ n ^7which is reached for K > 3. The peak is at K pea k = a/6/11 ~ 0.74, which 
corresponds roughly to the energy injection scale at K — kL/2ir — 1. The peak amplitude is 
smaller by a factor ~ 6 compared to the amplitude obtained when naively extrapolating the 
A; -11 / 3 behaviour down to the energy injection scale K — 1, as shown in Fig. |4| This is an 
important point when we compare the amplitude of the GW signal with previous estimates 



in the literature. Note that expression (28) is analytic for k — ► 0, which is required for 



causally generated, incompressible turbulence in the cosmological context [17]. 



2 In this paper we neglect the presence of a helical component in the velocity and magnetic field power 



spectra (c./. Eq. (37 1). Non-zero helicity, possibly arising from a macroscopic parity violation in the early 
universe, affects the decay of MHD turbulence, as demonstrated for example in [37], and the subsequent 
generation of GWs. For GW production by primordial helical MHD turbulence we refer to the analysis of 
Refs. [20,24]. 
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Figure 4: Comparison between the turbulent velocity power spectrum obtained by intersecting 
the k 2 behaviour at small scale with the inertial range A; -11 / 3 behaviour (blue), as done in the 
literature, with the Von Karman spectrum (black). 



In the inertial range, the characteristic velocity on a given scale 2ir/k is approximately 
given by [30] 

v\ ~ kE(k) ~ 6nC v ^K- 2/3 for K > 3 . (32) 



n 



rati 



Thus, the characteristic velocity on the energy injection scale, which is related to the total 
kinetic energy in the turbulence, is: 



v 2 T ~ 



6nC v 



On 



n 



rad 



AttC v (v 2 ). 



(33) 



3.3 The magnetic field power spectrum 

The power spectrum of the MHD processed magnetic field is closely related to the one of the 
turbulent velocity field. Here we treat the two sources exactly on the same footing. Since 
the energy density of a cosmological magnetic field scales like radiation, we can use Eq. ( 14 ) 
to evaluate the GW spectrum sourced by the magnetic field, provided that we define the 
normalized magnetic field vector 



(34) 



167rp ra d 

so that the transverse traceless (TT) part of the magnetic field energy momentum tensor is 



(TT) 



Bi{*,t)Bj{*,t) 
4tt 



(TT) 



4 



PmdOO [&i(x,t)6j(x,t)] 



(TT) 



(p+p)n*(x,t). 



(35) 

The normalized magnetic field bi is equivalent to the dimensionless turbulent velocity field 
i>j. We define a parameter analogous to Eq. (30), given by the ratio of the magnetic field 
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energy density to the radiation energy density 

Iu2\ 3 _ Pb ... (B 2 ) 

(b) = — = 2—— > Wlth = • 36 

2 i W P+P 87T 

In the following we assume equipartition between the magnetic and turbulent energy densities 
at the time when the latent heat is released, therefore (v 2 ) ~ (b 2 ) (however, the scaling of the 
GW spectra with the source energy density is kept explicit). The slope of the high frequency 
tail of the magnetic power spectrum in fully developed MHD turbulence is not precisely 
known: it could be of the Kolmogorov type, or it could satisfy the Iroshnikov-Kraichnan [58] 
or Goldreich-Sridhar [59] spectral slopes. While in the presence of a strong background 
magnetic field, the spectrum of the field component perpendicular to the background field 
is of the Iroshnikov-Kraichnan type according to Refs. [60-62], in the isotropic case, the 
relevant one in cosmology, the simulations of Ref. [61] indicate a Kolmogorov-type slope. 
However, to diversify the treatment of the magnetic source from the turbulent one, we 
choose to consider the Iroshnikov-Kraichnan spectrum; the GW spectrum resulting from a 
Kolmogorov magnetic field is not very different and can be readily derived from the turbulent 
one. 

The low frequency tail of the spectrum is determined by causality and by the fact that 
B is divergence free. Like for turbulence, we use the interpolating formula from Ref. [56] to 
find the equal time magnetic power spectrum 

(&i(M)&;(q,*)) = (2n) 3 6(k-q)P lJ P b (k,t), (37) 
PfU*\ 3 n ^b r3 K 2 / 1 for0<^<f ( . 

n(M) = 2 Cb n^ L oTIcW 1° for^>f. (38) 

The normalization constant C& = 16 J 3/2 ~ 0.0265 is calculated in the same way as in 



Eq. (31). The magnetic field and the turbulent flow being generated by the same physi- 
cal process, namely bubble nucleation and collision, we assume that they share the same 
correlation scale L and Kolmogorov microscale A. 

3.4 Freely decaying turbulence 

The turbulence stirring time is given by the duration of the phase transition which is typically 
much shorter than one Hubble time. Calling t in the time at which the phase transition 



starts, one has <C rfc~ n . As discussed in Section 3.1 , the typical time interval over which 
turbulence is established is the eddy turnover time tl- Moreover, Kolmogorov turbulence 
can be generated only if tl < rti n . This condition translates into a lower bound for the 
turbulent fluid velocity, in terms of the phase transition parameters: 

v f >(n- m /(3)v b . (39) 

If tl < T^in 1 ) turbulence sets in after a short interval of time At = tl < t- m . After the 
completion of the phase transition (once bubbles have percolated), the stirring is over and 
the turbulence enters the free decay regime, i.e. the total kinetic energy is dissipated (see 



for example [30]). The physical quantities appearing in the turbulent (28) and magnetic 
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(38) spectra are comoving, but they also have an additional 'absolute' time dependence due 
to the evolution of the turbulent cascade and the decay of the total energy. As already 
mentioned in the introduction, the existence of a maximal correlation length in the universe 
(the horizon) and the consequent Batchelor spectrum (i.e. k 2 at large scales), imply that 
(f2r/f2 ra( i)I/ 5 is constant in time [17,41,42]. According to this, we assume the following laws 
for the growth of the correlation scale and the decay of the kinetic energy: 

L(t) = L*(^— ^Y, 7 >0 (40) 




^rad 

When the phase transition starts, at i in , both the correlation length and the kinetic energy 
vanish. This insures that the source is continuous in time. At a time = t in + tl, after 
the completion of the phase transition, turbulence is fully developed. The stirring scale 
is given by L*, and 3^ is the total kinetic energy in the turbulent fluid, normalized to 
the radiation energy at time t*. We assume that the energy cascade responsible for the 
Kolmogorov spectrum starts at this stage. At times t > £ # = t in + t l turbulence enters 
the free decay phase, and the correlation scale and the kinetic energy evolve following the 
condition (fi-r/firad)-^ 5 =constant. 



The simple linear interpolation between these two behaviors given in (41) has been in- 
troduced to mimic the turning on of the source in the bubble collision case, inferred from 
numerical simulations, see [13]. The linear increase of the magnetic energy density in MHD 
turbulence has also been observed in simulations [44]. After this initial phase, for Kolmogorov 
turbulence which we shall adopt here, the energy decay law infers the value 7 = 2/7 (see 
e.g. [30] and the introduction). For generality, we keep 7 unspecified in the analytical 
formulae. When numerical results are presented, we substitute the value 7 = 2/7. 

In order to show explicitly its time evolution, we re-express the velocity power spec- 



trum Eq. (28) in terms of the time-independent variable K* = kL*/2iT. Introducing the 



dimensionless time variable 

t t; n 

y = , 

consequently L(t) = L^y 1 , and using K = K^y 1 we obtain: 

< y w if0<i,<i m d0<A-, <ij 

1 u 11 -a* ^ x{y) 



The time-dependence of P v (K^,y) is shown in Fig. [5j The initial phase, < y < 1 is 
inserted so that the source of GW increases smoothly from zero. As already discussed, 
the continuity of the source at initial time is an important issue for the resulting GW 
spectrum [14]. However, since the duration of this initial phase is short, we assume that 
the details of the source spectrum during this phase are not relevant and we do not model 
them in any detail. In particular, we do not expect the inertial range K* for K * > 3 
to be already developed in this initial phase, because the real energy cascade has not yet 
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started. Nevertheless, for simplicity we keep the same form of the spectrum as a function of 
K* in the two phases: the turbulent free decay phase, y > 1, is actually the most relevant 
one for the GW generation. During this phase, the large scale part of the power spectrum 
^ C 1 remains constant. We also show the characteristic velocity and the kinetic energy 
as functions of for different times in Fig. |6j 

In the following we assume that the turbulent magnetic field also undergoes the same 
time decay as the turbulent velocity field. Even though the precise decay law of the magnetic 
field energy density is not known in general for MHD turbulence, both analytical [32-35] 
and numerical [37-39] analyses seem to agree with the fact that the magnetic field power 
spectrum is persistent on large scales. Together with the condition that the spectrum, in 
a cosmological setting, should be of the Batchelor type (i.e. k 2 at large scales) due to the 
presence of a cosmological horizon, constancy in time at large scales entails the existence of 
a conserved quantity analogous to Loitsyansky's invariant: (Qb/^t^L 5 =constant. This 
gives the same decay as for the turbulent flow. Equivalent scaling laws (once generalized 
to the Batchelor case) are obtained from the argument of self-similarity [32,38] and direct 
cascade [34,39]. Therefore, we also assume in the magnetic case, 



if < y < 1 and < < ^ 
if y > 1 and < < ^ 
if > #r . 

Our argument for assuming this particular form of the time decay in MHD turbulence relies 
on approximative, analytical considerations, and given the high non-linearity of the problem 
it is conceivable that only numerical simulations will be able to find the correct scaling. For 
simplicity, and in order to be able to proceed with our analytical estimate, we are forced to 
make the above mentioned, rather crude assumptions for the scaling. Nonetheless, we would 
like to stress here once again the importance of the presence of a causal horizon, unavoidable 
in the cosmological setting. This prevents the formation of long range correlations at least 
beyond the horizon scale, a feature which certainly affects the decay law and that can not 
be accounted for in numerical MHD simulations which go on for times larger than the box 
size. 



n 
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rad 



(t) 
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^rad* 
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3.5 How long does turbulence last? 

In this section we confirm that turbulence is generated during a phase transition, and we 
determine when it ends according to the free decay picture described above. For this we eval- 



uate the Reynolds number defined in Eq. (23) at the energy injection scale L corresponding 



to K = kL/2ir = 1. We distinguish the physical length with a subscript p 

To ( go 
T \g{T) 



L P (T) = L(y) T 4 ( -^Y /3 (42) 
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K, = kLJ2n 



Figure 5: The normalized velocity power spectrum as a function of wavenumber K* for different 
times. Left: the phase in which the turbulence is developing, < y < 1. Right: the phase of free 



decay, y > 1. The Kolmogorov microscale is outside the plot range (c./. end of section 3.5) 




Figure 6: Left: the characteristic velocity v\ ~ kE(k) = 4irK^P v (K*) / Ll as a function of 
wavenumber at different times in the inertial range > 3 during the free decay phase, y > 1. 
Right: the kinetic energy E(k)/L* = 2K^P V (K*)/ Lf. The Kolmogorov microscale is outside the 
plot range. 



from the comoving scale L(y) = L^y 1 . The kinematic viscosity v{T) is derived in Ap- 
pendix |Bj 

( 22 T- 1 T > 100 GeV 

v{T) w I 5 10 s GeV 4 T' 5 T < 100 GeV (43) 
{ 2 10 9 GeV 4 T~ 5 T < 100 MeV 

The jumps in the viscosity introduce an uncertainty in the evaluation of the parameters, for 
which we can only give the correct order of magnitude. 



For vl = Vk(k = 2tt/L), we use the relation (32). Eq. (32) is valid only during the 
cascade, the free decay phase y > 1, and in the inertial range K > 3. To estimate the 
Reynolds number, we extrapolate it to K = 1, thus making a small error. However, the 
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following estimate is valid only for times y > 1. We have 



vl(y) ~ 67rC v ^y-^ , for y > 1 . (44) 

» 'rad* 



Inserting this in (23) we obtain for temperatures T < T* corresponding to y > 1 

Re(L(T)) = Re^)^^^ 1 ^^! 17 with ( 45 ) 



T z/(T) VsCO 
M^*) = Wfcrf^-^?? teY' 3 . (46) 



Re(L(T)) is plotted for different values of T* in Fig. [7j If we consider the EW phase 
transition, for which T* ~ 100 GeV, and we fix P/H* = 100, fir*/f2 ra d* = 2/9 corresponding 
to Vf = c s , Vb = 0.87 corresponding to detonations (see Sec. |3.1[ ), and 7 = 2/7, we find the 
value 

Re (L(T, = 100 GeV)) ~ 10 13 . (47) 

This confirms that turbulence develops once the primordial fluid is stirred on the scale 
= 2vb/P- In Appendix [B| we also derive 



-^-J , 1 MeV < T < 100 GeV, (48) 
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R m (L(T)) = Re(L(T))P m (T) ~ 10 12 Re (L(T)) f ^ ) . (19) 



so that 



V T 

R m (L(T, = 100 GeV)) ~ 10 17 . (50) 

We now want to determine the wavelength up to which the cascade is present, and the 
time up to which MHD turbulence persists. The turbulent cascade stops when viscosity 
becomes important, and this happens at scales smaller than the Kolmogorov microscale A, 
defined by 

Re (A) = ^ = l hence ^ = — Re (L) = Re (L) 3/4 , (51) 
v A v L 



where for the last equality we use (32). Therefore, the Kolmogorov microscale also grows 
during free decay, according to: 



T v(T) /sOZT 1/3X3 ' 



17. 



W = x -It,M{V ] 1 » TT - (52) 



The temperature dependence of A(T) and L(T) is illustrated in Fig. |8} The ratio of the 
initial scales determining the extension of the turbulent inertial range is L*/A* = O(10 10 ), 
for T* = 100 GeV, (3/H* = 100, fi T */^rad* = 2/9, v b = 0.87 and 7 = 2/7. 
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T (GeV) 

Figure 7: Evolution of the Reynolds number at the energy injection scale L P {T) as a function of 
temperature, assuming a phase transition with (3/l~t* = 100, = 2/9 and 7 = 2/7, for three 

phase transition temperatures: T* = 100, 10 6 and 10 8 GeV. The temperature at which the Reynolds 
number decays below 1 (horizontal line) represents approximatively the end of the turbulence: this 
typically happens at temperatures ~ 100 MeV, 2 GeV and 8 GeV respectively. 




Since A(T) grows faster than L(T), there always exists a temperature at which the two 
scale cross, as shown in Fig. [8] We define the end of turbulence when the entire inertial 
range (K > 3) is dissipated, namely when the dissipation scale has grown to reach 

^4=3 =► Re(L(T fin )) = 3 4 / 3 . (53) 

A(Jfi n J 
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If turbulence starts at T* = 100 GeV, using again the values f3/?i* = 100, QT*/^rad* = 2/9 
and 7 = 2/7, we find T^ n ~ 120 MeV. Therefore, turbulence acts as source of GWs for many 
Hubble times. 

Once turbulent motions are dissipated, the magnetic field fluctuations at scales larger 
than the dissipation scale remain frozen in the primordial plasma. Hence, in principle, the 
magnetic field continues acting as a source of GW also after the end of turbulence. However, 
we neglect this extra contribution to the GW spectra since it will not affect the spectrum 
at the interesting scales, around the peak, where the signal may be visible. These scales in 
fact have entered the horizon well before the end of turbulence and further GW production 
is strongly suppressed at later times. 



3.6 Time de-correlation of the spectrum of magneto-hydrodynamical 
turbulence 

In order to evaluate the GW power spectrum generated from MHD turbulence, the equal 



time velocity and magnetic field power spectra given in Eqs. (27, 28) and Eqs. (37, 38) are 



not enough. We need to know the unequal time velocity power spectrum, 
( Vi (k, *i)«;(q, h)) = (27r) 3 5(k - q) P v (k, h, t 2 ) , 



(54) 



and equivalently for the magnetic field. To model the unequal time velocity power spectrum, 



we multiply the equal time one (28) with the exponential time de-correlation proposed by 



Kraichnan in [43] and also used in [18]. The temporal de-correlation of the velocity field 
described in [43] operates only in the inertial range and during the cascade, which in our case 
means K > 3 and t > t*. The characteristic de-correlation time scale is given by the eddy 
turnover time: T£ = £/(2ve), for a characteristic eddy of size I = 2ir/k well in the inertial 
range. To model the de-correlation Kraichnan proposes a Gaussian functional form, 



giUM) = exp(-7r(t 1 - t 2 ) 2 /(4r/)) . 



(55) 



In our approach, the magnetic field undergoes the same de-correlation as the turbulent field. 

In the case of freely decaying MHD turbulence under consideration here, the eddy 
turnover time Ti is itself time-dependent, through the characteristic velocity on the same 



scale, ve. Using the formula of the characteristic velocity in the inertial range (32), and the 



evolution equations for the kinetic energy and the correlation scale (41) and (40) we obtain 



n{y) 



2v e 



t» y~ 



where 77? 



2\/6nC v 



(56) 



and the above equation is valid only during the cascade y > 1 and K > 3. Setting y — (t\ — 
tm)/T~L, and z = (t 2 — im)/T£, we arrive at the following expression for the velocity unequal 
time power spectrum (the magnetic field spectrum is readily derived from the expression 



21 



below, substituting C v by C b , Qt by fls and 17/6 by 11/4): 

p -<***> = ^^'''^'' fS^ (57) 



x 



for 2) < 3 

exp(-f(y-2:) 2 (^) 2 ) for 3<K(y, 2 ) < f and y, z > 1 



for K(y, z) > 



Here the notation ^^{y^z) (and so on) is to remind that the variables K, L, CIt/QtsA and 
T£ depend on time. In principle, we need to specify at which time these variables have to be 



evaluated. This choice must satisfy the constraint that the unequal time power spectrum ( 54 ) 
is symmetric under the exchange of t\ and t 2 - As explained in the next section, we avoid 
the problem by modeling de-correlation of the Kraichnan type directly in the anisotropic 
stresses of the source. 



4 The anisotropic stress power spectrum 

In order to determine the GW energy density power spectrum generated by the MHD tur- 



bulent source, we need to calculate the anisotropic stress power spectrum (10) at different 



times, substitute it into Eq. (14), and evaluate the time integral. 

The tensor anisotropic stress is the transverse traceless part of the energy momentum 
tensor E[y(k, t) = (PuPj m — l/2P ij -Pj m )Tj m (k, t), where we denote the dimensionless energy 
momentum tensor with a tilde, see Eq. ([6j. The part of the energy momentum tensor of our 
sources which contributes to the anisotropic stress is given by T) m (x, t) = f/(x, £)u m (x, t) 
and X^(x, t) = 6j(x, i)6 m (x, t). The anisotropic stress power spectrum is then 

(n^MOn*-^)) = (27r) 3 5(k- q)II(fc,ti,t 2 ) = V abcd {T ab {^M)T* c MM)) , (58) 
where we omit the tildes for simplicity and 

V abcd = \ PiaPjb - IPijPab^J (k) yP ic P jd - \PijPcdj (q) • (59) 
For the turbulent source, one has 

(n u (k 1 ,t 1 )n* J (k 2 ,t 2 )) = 

/d v f do 
(2tt)3 J (27r)3^ kl ~P' tl H(P'* 1 )^( k 2-q > *2M(q,*2)). (60) 

In order to proceed analytically, we assume that we can decompose the four point correla- 
tion function into products of the power spectra, using Wick's theorem like for a Gaussian 
random field. Since turbulence is not truly Gaussian, this is of course not strictly correct 
but it is usually adopted as a reasonable approximation to close the hierarchy (i.e. to avoid 
using equations involving higher order correlators). We refer to [21] and Appendix [C| for 
details of the determination of U(k,t 1 ,t 2 ). The bottom line is that the anisotropic stress 
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power spectrum is given by the convolution of the unequal time velocity power spectrum, 



c.f. Eq. (113). This quantity is in principle determined once we know how the source behaves 



in time. In our case it is given by the Kraichnan de-correlation model, once we have chosen 
an appropriate way of symmetrizing in time, see Eq. 



(57) 



On the other hand, the anisotropic stress power spectrum should, by definition, be a 
positive kernel, i.e. such that 



J dt x J dt 2 U{kMM)f{ti)r{t2)>U V/(t). 



(61) 



This follows simply from the definition of the power spectrum, Eq. (58), making use of the 
inequality 



/ d^iiyOMO/ftOf) >0. 



(62) 



Comparing Eq. (61) with the definition of the GW spectrum Eq. (14), we see that in the 



GW case the function f(t) is replaced by the Green function of the GW wave equation, 
f{ti) = cos(Mi)/ii or sm(kti)/ti (the factor l/t\ is absent for the short lasting case). Since 



the trigonometric functions are a complete basis, the property (61) is automatically satisfied 



if it holds for the Green function with arbitrary values of k. Therefore, inserting the unequal 



time source power spectrum Eq. (57), which accounts for Kraichnan de-correlation, in the 



convolution given by Eqs. (60) and (113), and performing the integration over the momenta, 
should give back an anisotropic stress power spectrum which is a positive kernel. However, 



we have tried different ways of symmetrizing Eq. (57) without succeeding in obtaining a 



positive kernel. This may be related to the inaccuracy of our analytical estimates, to our 
choice of symmetrization, or it could be an indication that either Wick's theorem or free decay 
together with Kraichnan de-correlation, although reasonable assumptions for the evolution 
of MHD turbulence, are not entirely appropriate for the source we are considering. Only a 
full numerical simulation of the turbulent and magnetic fields could help to understand the 
shortcomings of our approach. 

The same problem can arise when trying to obtain the source power spectrum by Fourier 
transforming the space correlation function, but this is circumvented since the source is 
modeled directly in Fourier space (the velocity power spectrum is well known in the MHD 
turbulent case). However, in order to calculate the GW spectrum one has to evaluate also 
the time Fourier transform: this is not a common procedure and no general solution is given 
in the literature. The problem is worsened by the fact that we want to account for the free 
decay of turbulence. If the only time dependence of the turbulent power spectrum was the 
Kraichnan de-correlation, which is Gaussian in the time difference, the spectrum resulting 
after time Fourier transforming would be positive. However, the free decay introduces an 
additional absolute time dependence which is more difficult to handle. Moreover, having to 
evaluate the anisotropic stress power spectrum further increases the difficulty of finding the 
correct time behaviour that would provide a positive kernel. Note that the same problem 
arose in the analytical evaluation of the GW signal coming from bubble collisions [12]. 

In order to proceed, we model the source in such a way, that H(k, t±, £2) is a positive kernel 
by construction. This is most easily done directly for the anisotropic stress power spectrum. 
In previous analyses [12] we have already tackled this problem and proposed three forms for 



the unequal time anisotropic stress power spectrum which are positive kernels. In Sec. 2.3 
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we presented the incoherent (17) and coherent (20) cases, for which the source is never 



correlated, respectively always correlated in time. We apply them in the following to the 
MHD turbulent source. We believe, however, that the top hat correlation introduced in 
Ref. [12] is the relevant one for MHD turbulence, since it best mimics the Kraichnan de- 
correlation. In the top hat model one assumes that Tl(k, t±, t 2 ) is correlated if \t\ — t 2 \ < x c /k 
and uncorrelated otherwise. Here x c is a parameter of order unity. We shall choose x c — 1 
for our numerical results. Since U(k,t 1 ,t 2 ) has to be symmetric in t\ and t 2 we set 



+ n(M2,*2)e(ti-* a )e(-^-(*i-* 2 ; 



(63) 



This way of correlating the source at unequal times is intermediate between the coherent and 
incoherent approximations; instead of being correlated at all times or only for t\ = t 2 , here 
we account for the fact that longer wavelengths de-correlate at larger time differences. While 
the term Q(ti — 1 2 ) is there only to make the function symmetric and does not influence the 
time continuity of the source, the term Q(x c /k — \t\ — t 2 \) should in principle be replaced 
by an exponential decay to keep the source continuous. We have tried this, and a part from 
a much slower convergence of the numerical integrals since the integrand oscillates rapidly 
for large values of k\t\ — t 2 \, we found no difference in the final result. Moreover, inserting 
(63) in the integral (11), we find that the GW energy power spectrum is not given by the 



Fourier transform of the source and continuity does not affect the final spectrum in this case 
(c.f. [14]). 

Assuming that the Kraichnan time de-correlation trivially extends from the turbulent 
velocity field to the anisotropic stress, it is clear that it gives a behaviour quite similar to 



the top hat de-correlation: according to (55), the source is no longer correlated for time 
intervals \t\ — 1 2 \ > (2/ y/n)re. Accounting for the fact that the eddy turnover time is simply 
the inverse of the characteristic frequency of the source ui = and that the GW Green 
function selects the diagonal of the time Fourier transform of the source, for which |u;| = k, 
we find that the Kraichnan de-correlation gives back the same condition as the top hat 
de-correlation, namely correlation is lost for time differences 



2 1 

TTUe 



X r 



with x c ~ 1 . 



(64) 



Even though the top hat case best reproduces turbulent de-correlation, in the following we 
evaluate the GW spectra also for the coherent and incoherent cases and compare the results. 
In order to proceed with the calculation, we now evaluate the equal time anisotropic stress 



power spectrum, see Eq. (63) 



4.1 The equal time anisotropic stress power spectrum for magneto- 
hydrodynamical turbulence 

The equal time anisotropic stress power spectrum is given by the convolution of the equal 
time velocity and magnetic field power spectra multiplied by an angular dependence coming 
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from the projector in Eq. (60), see Eq. (113). As already mentioned above, we refer to [21] 
and Appendix [O for a derivation. In terms of the variable K = K^y 1 ', using the velocity 



power spectrum given in Eq. (42), we find 



rati 



(y)) L 3 (y)l v (K,y,y) 



(65) 



where y = {t\ — tin) /tl and T v is given by 

Q 4 



Z v (K,y,y) 



n dQ \l + QY 7/e 



2 K 2 + Q\l + x 2 )-AKQ X 
+ K 2 — 2 X KQ + Q 2 ) 17 /6 



(66) 



with Q = Lq/2ii, x = k-q. The magnetic field anisotropic stress power spectrum is equivalent 
to the above expressions, substituting Q T with Q B , C v with C b and the power law 17/6 with 
11/4. The integrals cannot be performed analytically, so we solve them numerically and 
derive fits in terms of the variable K = K^y 1 . Within this approach we do not account for 
the small scale cutoff L/X when fitting the convolution. We explain below how this cutoff is 
taken into account. We find, for the turbulence and the magnetic field respectively: 



l v (K*,y,y) 
I b (K # ,y,y) 



0.098 



0.12 



4/3 



+ 



4/3 



3.3 
3.5 



11/3' 



7/2' 



-1 



(67) 
(68) 



The fits are shown in Fig. [9] 



4.2 The final time of the magneto- hydrodynamical turbulent source 



As demonstrated in Sec. |3.5[ MHD turbulence can act as a source of GWs for many Hubble 
times, while it undergoes free decay. The source generally switches off at the end of the 



turbulence, when Re(L p (T fin )) = 3 4 / 3 (see Eq. (53)). If T* = 100 GeV, we have found 

that the Kolmogorov microscale grows in 



fin 



120 MeV. However, we have seen in Sec. 3.5 



time, and for wavenumbers above this upper cutoff the source has decayed. This means that 



the final time of integration in Eq. (11) must be defined as the minimum of t(T nn ) and the 



time at which a given mode k is equal to the upper cutoff. For later times, that mode k is 
not generating GWs any longer; however, if the time at which k is equal to the upper cutoff 
comes after the end of the turbulence, we should take the end of the turbulence t(T fin ) as 
the final time of action of the source on the scale k. This is the way in which we include the 



upper cutoff L/X appearing in Eq. (57) in the calculation of the GW spectrum, even though 



we have neglected it for simplicity in the evaluation of the anisotropic stress power spectra 



(67), (68). 



We introduce % as the time at which k = 4ir/X(tk), where the extra factor of 2 comes 



from the fact that I v (K*,y } z), being the convolution of the velocity power spectrum (57) 
(and equivalently for the magnetic field power spectrum), it will go to zero at twice the 
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Figure 9: Black, solid lines: the integral in Eq. (66), left: for the turbulence anisotropic stress, 
and right: for the magnetic field anisotropic stress, as a function of K = K^y^, together with their 
fits given in Eqs. (67) and (68) (shown as red, dashed lines). 



velocity power spectrum cutoff. Using the time evolution of the dissipation scale A, Eq. (52), 
we find 



2 f t>^ \ 



A* K* \t l J 



- 



177/8+3 



if T* < 100 GeV , T fin > 100 MeV 



(69) 



A* K* 



8 
17- 



if T fin > 100 GeV , (70) 

where the different behaviour depending on the initial and final temperature is due to the 
different evolution of the viscosity, given in Eq. (43). Finally, the time at which turbulence 
at a given scale k ends is tsu(k) given by 



Un(k) 



mm 



tl 



fin J 



tl 



tl 



(71) 



where t(T fin ) — T / (T fin H ^Q rSLd ) and T fin denotes the end of the MHD turbulence by the 



dissipation of the entire Kolmogorov range, denned in Section [375] For the usual set of values 
T* = 100 GeV, /3/H* = 100, fi T */ft rad * = 2 /9, v b = 0.87 and 7 = 2/7, we find 



*fin(fc) 



2 x 10 4 



tl 



28 

2 10 1 4 \ 101 



for K* < 0.07 
for K* > 0.07. 



(72) 
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5 The gravitational wave spectrum 



We are now ready to evaluate the integrals in Eq. (Jill). We consider the three approximations 
for the unequal time correlator of the anisotropic stress tensor mentioned above, namely 
incoherent, coherent and top hat. All we need is the anisotropic stress power spectrum 
taken at equal times, which is given in Eqs. (67) and (68). These are then inserted into 



Eq. (11). The GW power spectra obtained in this way are always positive. The figures of 



the GW energy density spectra shown below are calculated for the set of parameter values 
given at the end of the previous section. Under the equipartition hypothesis, we fix the 
magnetic field energy density to the same value of the turbulent one: -^ XjL - = = | . 



Even though we have demonstrated in Sec. 3.5 that MHD turbulence can last for many 
Hubble times, at the end of this section we also consider the case of GW generated by MHD 
turbulence confined in time to the duration of the phase transition. We find this analysis 
illuminating to understand our results and useful to compare with the results obtained in 
previous works. 

From now on we use a common notation for both the GW spectra generated by turbulence 
and by the magnetic field, with s = v or b denoting respectively the turbulence and magnetic 
field source. 

• Incoherent approximation: in this case the anisotropic stress spectrum at unequal 



times is given by (c.f. Eq. (17) and (67) resp (68)): 

n(K*,y,z) = \*Cl (j^(f)) L 3 (y)I s (K*,y,y) 8{y - z) 



(73) 



where we have chosen tl as the short characteristic time over which the source remains 
coherent. From the anisotropic stress formula given above we find the spectrum: 
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(74) 



l s (K*,y,y) 



'o [y+%] ' ' ' h [y + %] 
where tl = L*/(2vl) is the initial eddy turnover time at the scale for y = 1, and 



v L 



i'tnCvTT^f- denotes the initial eddy turn-over speed, Eq. (|44|); moreover, see Sec. 

1 where 



3.4 



tin 



tl 



v_lI_ 



(75) 



?/ nn is given by (tf\ n (k) —t- m )/rL and t& n (k) is defined in Section 4.2 The spectra X s (K # ,y,y) 



are given in Eqs. (\67v and (68). 



We recover the same behaviour as in section 2.3 The GW power spectrum is proportional 
to the phase space volume K% times the source power spectrum. The slope at large scales 
is therefore the one at small scales is K 3+n where n = —11/3 for the turbulence and 
n = —7/2 for the magnetic field. The peak is at f^P cak ~ 5.9 for turbulence and Kf eak ~ 7 
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for the magnetic field. The results are shown in Fig. 11 for the usual choice of the parameter 
values. The incoherent approximation is the one leading to the highest peak amplitude: 



dQawhg 
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dlogk 
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(mag) 


dlogk 


^-peak 



~ 4 x icr 10 

~ 3 x 10~ 10 



Kl eak -5.9 

^-peak ^ j 



(76) 
(77) 



• Coherent approximation: According to Eq. (20), we have 

Il{K*,y,z) = |ttC 2 §^(y)Ll(y) ^( z )Li{z) y/l s (K„y,y)y/l s (K„z,z) 

^ » 'rad " 'rad 

The GW spectrum is positive by construction, but it is oscillatory: 



(78) 
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dy V t . y/l s {K*,y,y) cos 
y+ir V V L 
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37/2+1 



7 7 /2 



y/l a (K*,y,y) sin ( 



*2/fin 



c?y ^ + ^ V x s( K *,y,y) sin 



As discussed in Ref. [14] and in Section 2.3, the calculation in the coherent case is more 
involved than the incoherent one, since the GW power spectrum is not simply proportional 
to the source power spectrum, but it is given by the square of its time Fourier transform. 
The source is characterized by the space correlation scale and the time correlation scale 
tl, related by = 2vlTl- Since vl < 0.4 (the upper bound is given by Eq. (33) with (v 2 ) = 



1/3), we have < T£. On scales larger than both the characteristic spatial correlation scale 
and time correlation scale tl, the Fourier transform of the source is constant because the 
source is not correlated (white noise). Therefore, for wave-numbers k < 2tt/tl < 2ir/L*, 
corresponding to K* < L*/tl, we recover the K% behaviour like in the incoherent case. 
However, for k > 2tt/tl the time Fourier transform is no longer constant and starts to 
decay as a power law, the exponent depending on the time differentiability properties of 
the source [14]. Since we have chosen a source which is continuous but not differentiable at 
initial time i m (see the time evolution laws of the energy and the correlation scales given 



in Eqs. (40, 41)), this implies a decay like k 2 for the time Fourier transform of the source. 



Therefore, the GW power spectrum (79) (the square of the Fourier transform, multiplied by 
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Kl) decays like K~ x at intermediate scales L*/tl < K* < 1. This behaviour is satisfied up 
to the wave-number corresponding to the spatial correlation scale: k ~ 2ir/L*. Afterwards, 
the spectrum decays with a power law given by the time Fourier transform multiplied by the 
power law decay of the source. In the case of turbulence, for K* > 1 one has therefore the 
power law decay K* 14 ^ 3 = K% x (K~ 2 ) 2 x (K* 11 ^ 6 ) 2 ; while in the magnetic field case this 
becomes K~ 9/2 = Kl x (K~ 2 ) 2 x (K; 7/A ) 2 . 

The form of the power spectrum, in particular the wave-number at which the spectrum 

we show the GW spectrum 



o r : • in Fig. 



10 



peaks is determined by the ratio 
from turbulence for the coherent case for two values of this ratio: the case when the kinetic 
energy density in the turbulence in maximal (v 2 ) = 1/3, corresponding to -^j- = 2/9, and 
another one in which the kinetic energy involved is much smaller, (v 2 ) — 1 n_ ' ! 



to fer = 2/3 x 10- 



10 6 , corresponding 
3 . Notice that not only the amplitude but also the peak position differs, 



and that the power law behaviour derived above is recovered. The form of the GW spectrum 
in the coherent case is also shown in Fig. 11 for the usual choice of the parameters. The 
coherent approximation leads to the smallest peak amplitude: 
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Figure 10: The GW energy density power spectrum from turbulence in the coherent case for 
% = 100 GeV, P/H* = 100, v b = 0.87, 7 = 2/7, and 5fm = 47.75. Blue, solid: T */Q rad * = 2/9. 
The peak is at about K^ ea ~ 0.4, not far from the value L*/tl = 2{Q-kC v ^It* /^^d*) 1 ^ 2 — 0.8. The 
small scale behaviour for K* > K* cak is 14 ^ 3 . Red, dashed: /^rad* = 2/3 x 10~ 3 . The peak 
position is at ~ 0.02, again not far from the value L^/tl = 2(6itC v £It* I '^rad*) 3 — 0.04. In this 
case the slope K~ l for intermediate wave- numbers is well visible. 
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Figure 11: The GW energy density spectrum in the incoherent (red, long-dashed), top hat (black, 
short-dashed) and coherent (blue solid) cases. Left from turbulence, right from magnetic field, for 
T* = 100 GeV, 3/H* = 100, = 2/9, v b = 0.87, 7 = 2/7, g &n = 47.75 (and x c = 1 for the 

tophat case). 



• Top hat approximation: this is the most realistic case for the MHD turbulent source, 
since it mimics a de-correlation in time of the Kraichnan type as discussed in Section |4| 
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We choose the value x c = 1, as in the Kraichnan model, so that the integral determining the 
GW spectrum is positive. It is given by: 
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(84) 



As in the incoherent case, the spectrum bears no relation with the time Fourier transform of 
the source (c.f. [14]). The integral in z can be estimated simply as the integrand evaluated 
at the lower bound, multiplied by one oscillation period: therefore, at high wave-numbers 
we expect the slope K* 5 ^ 3 for the turbulent case, and K* 3 ^ 2 for the magnetic field case. 
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Moreover, the peak position corresponds to the spatial correlation scale of the source: ~ 
3.5 for the turbulence and K* ~ 3.7 for the magnetic field. The result is shown in Fig. [IT] and 
12, for x c = 1; the dependence on the value of < x c < ir is weak. In this case the amplitude 
at the peak takes an intermediate value between the incoherent and coherent cases: 
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Figure 12: The GW energy density spectrum in the top hat case. Blue, solid from turbulence 
and red, dashed from the magnetic field for % = 100 GeV, P/H* = 100, jpa*- = 2/9, v b = 0.87, 
7 = 2/7, g &n = 47.75, x c = 1. 



6 Discussion 

6.1 Comparison with the gravitational wave spectrum from bub- 
ble collisions 

In the case of bubble collisions [12], we have also studied the different assumptions for the 
unequal time power spectrum of the source: coherent, incoherent and top-hat. Numerical 
simulations of bubble collisions [13] indicate that the coherent case is the relevant one for 
bubbles. This can be understood by the following argument: we consider a bubble collision 
event starting at time t n with tensor anisotropic stress given by / n (k, t — t n ). Summing over 
all the collision events of the phase transition we obtain for the total anisotropic stress of 
bubble collisions (c.f. [14]) 

N N 

(n(k, *)n*(k', o> = Yl E (e l(k - x "- k '- Xm) /n(k, t - t n )t(k>, t> - O) , (87) 

n=l m=l 
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where x n is the center of the n-th collision event. To recover the simulation result, we now 
have to make two fundamental assumptions. First we assume that different collision events 
are not correlated, 

(e i(k.x„-k'.x m )) = V -i 5nm 5(yi - k') , (88) 
here the volume V restores the dimensions. This leads to 

i N 

n(fc, t, = j— £</„(k, t - t n )f n (k, t' - Q) . (89) 

^ ' n=l 

Second, we assume that one typical collision event is totally coherent, so that (/ n (k, At)/*(k, At')) 
/(k, At)/*(k, At') where /(k, At) = a/ (|/n(k, At) | 2 ) is the square root of the power spec- 
trum of a typical collision event, and At = t — t n . Furthermore, since bubbles only exist 
during the phase transition the duration of which is much shorter than one Hubble time, the 
expansion of the universe can be neglected and the anisotropic stress spectrum is a function 



of the time difference t — t' = At — At' only. Equation ( 89 ) becomes 



2N 

n(fc, t, = TCTy /( k > At )/( k > At ') = Vn(fc,t,t)n(fc,t',f) , (90) 

which is the form of the coherent approximation. Therefore, assuming that one collision 
event is totally coherent in time and that different collisions are uncorrelated implies that 
bubble collisions represent a totally coherent source of GWs. 

However, for turbulence, we do not expect this to hold. There are no well isolated 
uncorrelated events which can be treated independently. Besides, one expects correlations 
to decay in time over a timescale which is related to the spatial extension of the source. 



This has motivated the Kraichnan de-correlation ansatz given in [43] and Eq. (55), which 



we have simplified to the the top hat de- correlation in Section [4] in order to obtain a positive 



kernel, see Eq. (63) and the discussion following it. To summarize, we can conclude that 
bubble collisions are well represented by the coherent case, while MHD turbulence is well 
represented by the top-hat case. 

6.2 Comparison with short-lasting turbulence 

Once MHD turbulence is generated, it decays following the 'absolute' time dependence de- 



scribed in Section |3.2[ which applies to freely decaying, non-helical turbulence. In this 
section we compare our results with the ones from MHD turbulence that lasts for less than 
one Hubble time and where the time-dependent decay is neglected. This was always assumed 
in previous analyses. This allows us to verify the general statements of Section 2J3 in a more 
realistic case. Note however that, with respect to previous analyses which assumed either a 
discontinuous [17] or a stationary [15,18,20] source, here we consider a source with a finite 
time duration, but continuous in time. The importance of having a continuous source has 
been discussed in [14]. 



The aim of this section is simply to test the results of section [273] in a realistic case. We 
therefore concentrate only on the source from the turbulent velocity field, and we do not 
discuss the magnetic field for which the results are similar. We set the duration of the source 
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to one eddy turnover time tl- t& n = tin + tl [15]. We model the switching on and off of the 
source with the same function as in bubble collisions [14]: 



f(y) = 4y{l-y) 



(91) 



with y = (t — t in )/TL. Since L and flr/^rad do not evolve in time, the turbulent spectrum 
becomes 
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where I V (K*) is given by Eq. (67) with y = 1. Inserting this in Eq. (13) one finds after some 
manipulations 
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Figure 13: The GW energy density spectrum for short lasting turbulence. Blue, solid: coherent; 
black, short-dashed: top hat; red, long-dashed: incoherent. 
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Figure 14: Comparison between the long-lasting case, represented by thick lines and the short- 
lasting one, represented by thin lines (incoherent: red, long-dashed; coherent: blue, solid; top hat: 
black, short-dashed). 



Fig. 13 shows the GW energy density spectra from short lasting turbulence in the inco- 
herent, coherent and top hat cases. The relative amplitudes between the different approxi- 
mations is similar to the long lasting case and the peak positions are also not very different. 
In Fig. 14 we compare the long-lasting and the short-lasting cases. 

One first notices that the change in amplitude of the GW spectrum is generically much 



less than expected from our general arguments in Section 2.3 This is because our realistic 



source, unlike the one considered in Section 2.3, is decaying with the decay time tl, which 
here corresponds to the duration of the short lasting source. Therefore, the fact that the 
source is long lasting does not amplify the signal as expected, because its characteristic decay 
time is the same as the duration of the short lasting source. In the incoherent approximation, 
the overall amplification of the long lasting source is about a factor two, instead of the factor 
{tlH*)~ 1 which we found for the toy model. The long lasting coherent approximation, on 
the other hand, is amplified at large scales by about a factor of {tlH*)~ 1 — 80 compared 



to the short lasting one, instead of the expected {tlH*) 2 ■ As already seen in Section 2.3 



in the coherent case the signal at the peak from short-lasting turbulence is higher than 
the long-lasting result: this comes from the fact that in the long-lasting case interference 
can reduce the final amplitude if the source is present and coherent over several oscillation 
periods. 

The top-hat case, which we did not analyze in Section 2^3 shows an intermediate behavior 
between the incoherent and coherent ones: at long wavelengths, <C K^ eak , the long lasting 
signal is enhanced by about two orders of magnitude, i.e. the same amplification as for the 
coherent case, (tlH*) -1 ; around the peak, the long and short lasting top hat cases differ by 
about a factor two, equivalent to the incoherent case. 

To summarize, due to the fact that the long lasting source decays with a characteristic 
time tl corresponding to the duration of the short lasting source, the amplification factor 
(tl7Y*)~ 2 , expected at large scales, is reduced to (t^W*) -1 for the coherent and top-hat 
approximations, and the amplification is virtually absent in the incoherent approximation. 
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In all cases it is a factor (tlTC*) less than expected due to the rapid decay of the source. 

7 Conclusion 

In this work, we have calculated the GW emission from MHD turbulence generated during a 
first order phase transition and freely decaying afterwards. This is the first paper which takes 
into account the free decay of turbulence, and models the source in a continuous fashion. 

For the source power spectrum we use a new ansatz that interpolates analytically between 
the large scale and small scale behaviors, determined respectively by causality and by the 
Kolmogorov (or Iroshnikov-Kraichnan) theory. Previous analyses had either considered only 
the Kolmogorov range, ignoring the large scale part of the spectrum and continuing the 
Kolmogorov slope up to the peak [15,16,18,20], or joined the two behaviors at the peak [17]. 
This caused an overestimation of the source spectrum amplitude at the peak of about a factor 
six, leading to an overestimation of nearly two orders of magnitude for the GW spectrum 
(c.f. for example Eq. (56) of [17]). The interpolating formula which we adopt here models 
the spectrum of MHD turbulence in a more realistic way. 

We also take into account the time de-correlation of MHD turbulence following the model 
proposed by Kraichnan in [43] . However, in order to recover a positive result for the GW en- 
ergy density spectrum, we cannot directly apply the Kraichnan de-correlation in the velocity 
and magnetic field power spectra, but we have to model it in Fourier space as a de-correlation 
of the anisotropic stress power spectrum. We claim that in the case of MHD turbulence, 
neither the coherent [17] nor the stationary [15, 18, 20] approximations previously used in 
the literature are the correct ones, but the anisotropic stresses at different times have to be 
modeled in a way similar to our top hat ansatz. 

Moreover, previous analyses have considered either a MHD source which is discontinuous 
in time (i.e. instantaneous turning on of the Kolmogorov spectrum) [17] or a source which is 
stationary (i.e. neglecting the fact that the source is actually turned on and off) [15, 18,20]. 
Here instead the source is continuous in time, starting with zero energy and building up 
the Kolmogorov spectrum after one eddy turnover time. It then starts free decay, since the 
stirring due to the phase transition lasts only for about one eddy turnover time. Due to 
the free decay, the source is absent on scales smaller than the time dependent Kolmogorov 
microscale, and is completely dissipated once the Kolmogorov microscale has reached the 
energy injection scale and the entire Kolmogorov range has decayed. We have evaluated the 
temperature of the universe, Tg n , at which this happens, as a function of the temperature at 
which MHD turbulence is generated, T*, and we have found that the process of dissipation 
lasts for many Hubble times. For instance, for the EW phase transition, T* ~ 100 GeV 
and T fin ~ 120 MeV. Therefore, MHD turbulence has to be modeled as a long lasting GW 
source. Nevertheless, since the characteristic decay time of the source is still given by the 
eddy turnover time which is much smaller than the expansion time of the Universe, the 
amplification due to the long duration is less significant than what is expected from a source 
which is not decaying in time. Especially the peak amplitude is enhanced only by about a 
factor of two due to the long duration of the source. On the other hand, the long duration 
of the source becomes important on very large scales: here the GW spectrum from the long 
lasting source is amplified with respect to the short lasting case by about two orders of 
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magnitude. 

The top hat ansatz together with the time continuity of the source have some interesting 
consequences on the peak position and amplitude of the GW spectrum. In Ref. [14] it has 
been shown that, in the coherent case, time continuity affects the slope of the spectrum at 
small scales and also moves the peak of the GW spectrum from the characteristic length scale 
of the source (here L*) to its characteristic time scale (here tl). If the two scales are well 
separated (for instance if the fluid velocity is significantly smaller than the speed of light), 
this causes a reduction of the amplitude at the peak. The fact that for MHD turbulence 
the top hat ansatz is the relevant one, fixes the peak position of the GW spectrum at the 
characteristic length scale of the turbulent source. 

To summarize, our final result for the GW spectrum from MHD turbulence for the most 
realistic case, the top-hat de-correlation, has the following main features: 

• The peak frequency is given by ~ 3.5 for the turbulence and K* ~ 3.7 for the 
magnetic field. Given the degree of precision of our analytical estimate, we can neglect 
this small difference. The peak frequency from MHD turbulence thus corresponds to 

f MHD _ fc pcak 3-5 _ J3^ 

/peak " 2n ~ ^ ~ ^ ^ W 
o in 2 tt ( 9* V /6 pi 

~ 3x KT 2 mHz — 

V100/ lOOGeVK^fe 

~ 3.4 mHz , 

where we have used H* = (go/g*) 1 ^ 6 (T*/Tq) ifov^rad, and the last line gives the value 
for the EW phase transition with 0/H* = 100, T* = 100 GeV, v f = l/y/3 and 
v b ~ 0.87. 

• The peak amplitude for the above values of the parameters is h^flGw{.K^ cak ) ~ 10~ n . 
Using that h(f) = 1.26 x 10~ 15 a/ hlQ GW (f)(mKz/ f) (see [63]), we obtain the maximal 



GW amplitude h(f peak ) ~2x 10~ 21 which is detectable with LISA (see Fig. 15). The 
suppression of the peak amplitude by more than one order of magnitude compared 
e.g. to Ref. [18], is mainly due to the fact that we use the more realistic interpolating 
spectrum for the MHD turbulent source (see Fig. [4]). 

The slope of dflcw I 'dlog(k) is the usual k 3 on large scales k < fc pea k, but on small 
scales, k > A; pC ak, it decays only like k~ 5 ^ 3 and k~ 3 ^ 2 for turbulence and magnetic field 
respectively, see Fig. [12} 



In Figs. 15 and 16 we compare the GW spectrum from MHD turbulence with the ex- 
perimental sensitivities of LISA [25], AGIS [64], LIGO [65] and the Big Bang Observer 
(BBO) [66]. Note that the sensitivity curve of LISA plotted in the figures represents only 
the noise of the instrument, and the data analysis can actually improve the detection down 
to a level of h 2 Q VL GW ~ 5 x 10~ 13 [67]. 

With respect to the GW signal from bubble collisions analyzed in Ref. [12], the peak 
frequency of the GW spectrum from MHD turbulence is larger by about a factor two. Note 
however that, as already pointed out in [14], the work of Ref. [12] has to be corrected in 
two aspects: first, the source analyzed there was not continuous, and secondly, the relevant 
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approximation for the anisotropic stresses generated by bubble collisions is the coherent and 
not the top hat one. These modifications are work in progress and may well lead to some 
correction in the results of [12]. The peak amplitude of the GW signal from MHD turbulence 
is somewhat higher than the signal from bubble collisions, but it decays faster than the k" 1 
decay which has been seen in the latest simulations of bubble collisions [13]. 

Nevertheless, there is still considerable uncertainty in our analytical modeling which 
probably can only be addressed by numerical simulations of relativistic MHD turbulence of 
the kind developed after a first order phase transition. For example, in this work we have 
added the results from the turbulent velocity field and the magnetic field incoherently. One 
could argue, however, that these fields are correlated and have both a Kolmogorov spectrum 
in the inertial range. One could then use Sj = hi + v i ~ 2v i as transverse vector-field with 
Iljj(s) ~ 4Iljj(v). Instead of the results presented here we would then obtain 16 times the 
G W energy density spectrum from turbulence, which would enhance the total result for flow 
by about a factor of 8 and the one for h(f) by \^8. Therefore, our results are probably to 
be taken within about a factor of a few for the GW amplitude h(f) and within an order of 
magnitude for the energy density dflcw I 'dlog(k) . 



T, =100 GeV, /?/H=100 T, =100 GeV, /?/H=100 




f (Hz) f (Hz) 

Figure 15: The GW energy density (left) and the GW characteristic amplitude (right) from 
turbulence only (blue, solid), magnetic field only (red, dashed) and the total MHD turbulence 
(black, solid) generated at the EW phase transition with T* = 100 GeV, (3/Tt* = 100, Qs* /^rad* = 
2/9, Vb = 0.87, 7 = 2/7, x c = 1 and ga n = 47.75, together with the sensitivities of LISA [25], 
BBO [66], AGIS [64] and 'BBO Corr' (improved from BBO with data analysis) taken from [68]. 
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Figure 16: Sensitivities of LISA, AGIS, BBO and Advanced LIGO (orange) compared with two 
GW spectra (black) generated by MHD turbulence from a phase transition at respectively T* = 100 
GeV with 0/H* = 100, and T* = 5.10 6 GeV with (3/H* = 50; Qs*/^«*i* = 2 /9, v b = 0.87, 7 = 2/7, 
and x c = 1. The Advanced LIGO sensitivity is optimized by making use of correlations between 
two ground-based detectors [69] . 



A Analytical expressions for Section 2.3 



Here we give the full expression for Eqs. (19) and (22) 
• Incoherent constant source 



F(t in , tftn, Af) 



(— ) 

Vflfin I 



<fln 
At 



( 90 \ 3 At 
V Sfin / *in 



I \jj 3 V ti n / 



,t fln -At/2 

long-lasting, 
short-lasting. 



(97) 
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Coherent constant source 



(— ) 

V<7fi„ / 



,3j_ 
> Sfin > Ax 2 



x fin Ci(x fin ) + x in Ci(xi 



F(%iny -^fin; Ax) < 



-(sfin - Ax/2)Ci(x fin - Ax/2) - (x in + Ax/2)Ci(x in + Ax/2) 
— sin(xfi n ) — sin(xi n ) + sin(xfi n — Ax/2) + sin(xi n + Ax/2) 

- x fin Si(x fin ) + x in Si(x in ) 
x fin - Ax/2)Si(x fin - Ax/2) - (x in + Ax/2)Si(x in + Ax/2) 

1 2 

cos(xfi n ) + cos(xi n ) — cos(xfi n — Ax/2) — cos(xi n + Ax/2) 



- (ft)" [(Ci(sfin) - Ci(x in )) 2 + (Si(x fin ) - Si(x in )) 2 ] + 0{Ax) 
long-lasting , 

fgo\\ 64(2vr) 2 sin 4 ((x fln -x' in )/4) 
V9* > xL 



(^fin ^in)^ 



short-lasting. 



B Viscosity and magnetic diffusivity 



The Reynolds number defined in Eq. ( 23 ) is inversely proportional to the kinematic viscosity 
v given by 

" = (98) 
P + P 

where fj is the shear viscosity. The kinematic viscosity is the transport coefficient that 
characterizes the diffusion of transverse momentum due to collisions in a medium, and is 
roughly the mean free path £ mfp of excitations. In fact one has [70] 



4 7T 2 



-mfp 

~5~ 



(99) 



The largest viscosity comes from the weakest interactions, since it is inversely proportional 
to the scattering cross section of the processes responsible for transport. Simple parametric 
estimates using kinetic theory show that the shear viscosity at high temperature (where T 
is much larger than the mass of the diffusing particle) behaves as (to leading-log accuracy): 

T 3 

V = c— - (100) 

where g is the appropriate coupling constant (depending on the temperature and the length 
scale at which one wants to compute the Reynolds number) and C is a numerical coefficient 
that can only be obtained from a detailed analysis. 

After EW symmetry breaking, neutrino interactions are suppressed by a factor (T/Mjy) 4 . 
In this regime, neutrinos have the longest mean free path and dominate the viscosity. We 
use [71] 

4nf P « (3GIT 5 )- 1 , (101) 
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leading to 



u(T < 100 GeV) « 4.9 x 10 f 



GeV 4 

J^5 



(102) 



At temperatures smaller than 100 MeV, after the QCD phase transition, the particle content 
changes and consequently the neutrino mean free path increases to 



10 



2 m5\-l 



leading to 



imfp » ~^( G f t 



v(T < 100 MeV) » 1.6 x 10 £ 



GeV 
y5 



(103) 



(104) 



At temperatures above the EW phase transition, neutrino interactions are no longer sup- 
pressed. The shear viscosity is dominated by right handed lepton transport and given by [72] 



5\ 3 2 /12\ 5 3/2 T 3 

V ~ [ 2j CU \n) 9tt 2 + 224(5 + 1/2) g'* log g'- 1 



(105) 



where g' is the hypercharge coupling. This leads to 



v{T > 100 GeV) 



21.6 



(106) 



The evolution of v with temperature is plotted in Fig, 17 The neutrinos remain the relevant 



kinematic viscosity 




T (GeV) 



Figure 17: Evolution with temperature of the kinematic viscosity v(T). 



particles controlling the viscosity until they decouple at T ~ 1.4 MeV, after which photons 
take over. Even if there was some source of turbulence after the EW phase transition, 
turbulence is expected to terminate anyway around 1 MeV. Indeed, e + e~ annihilation reduces 
the plasma electron population and increases the photon diffusion length hence also the 
kinematic viscosity, leading to a decrease of the Reynolds number below one. 
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Since we have only found non-relativistic derivations in the literature [31,55,73], let us 
estimate here in some detail the magnetic diffusivity and the magnetic Prandl number for 
relativistic electrons in the cosmic plasma with temperatures 1 MeV < T < 100 GeV, 

= R m (L, T) KT) 
ml J "Re(L,T) fi(T) ' [Wi) 

We want to determine P m (^) when the electrons are relativistic and their dominant inter- 
actions are electromagnetic. 

To determine the magnetic diffusivity yu(T) we derive an expression for the conductivity 
cr(T). The Lorentz force acting on an electron is 

m e — = eF lLV u v . 
dr 

If we average this equation over a fluid element containing many electrons, the magnetic field 
term is sub-dominant. Even though the electrons are highly relativistic, the average fluid 
velocity is small. Furthermore 7 = 1/ a/1 — fg — T ' /m e is nearly constant and we may neglect 
the contribution d'j/dr from du l /dr = d^v^/dr above. With dr = / y~ 1 dt = (m e /T)dt, this 
yields the following equation for the mean velocity of the electron fluid: 

dv e „ 
— = -E. 
dt T 

If we denote the collision time for the electrons by t c , they can acquire velocities of the order 
v ~ j, Et c between successive collisions. Hence the current is 

e 2 n e _, 

.1 ~ en P v ~ t r — — E = crE 



e 2 n P 



so that the conductivity becomes 

o = t. 

We now derive an estimate for t c from Coulomb interactions. For a strong collision between 
the electron and another charged particle we need an impact parameter b such that e 2 /b > 
E e ~ T. Hence the cross section becomes o t ~ nb 2 ~ 7re 4 /T 2 (this simple argumentation 
neglects the Coulomb logarithms which enhance the cross section by ln(l/a min ) where a m \ n 
is the minimal deflection angle [73]). With v e = 1 the time between collisions is therefore 
t c = l/{a t n e ) ~ T 2 /(vre 4 n e ) and 

«T~-^r. (108) 

ire z 

Note that this result is independent of the electron density. This is physically sensible as n e 
enhances the current on the one hand but it reduces in the same way the collision time. 



With ( 108 ) we obtain for the magnetic diffusivity 



1 e 2 10~ 3 , . 

PV J Ana AT T K J 
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Inserting the kinematic viscosity from Eqs. (102) or (104) we obtain for the Prandl number 



12 /GeVV , . 

^10 12 — . (110) 



fi { T ) 



This number is larger than 1 for all temperatures 1 MeV< T < 100 GeV where the derivation 
applies. Hence currents and magnetic fields can develop and we are in the regime where MHD 
turbulence applies. 

C The unequal time anisotropic stress power spectrum 



As mentioned in the main text, the unequal time anisotropic stress power spectrum H(k, ti, t%) 
is given by the convolution of the unequal time source power spectrum. The details of the 
derivation (for the formally identical case of a magnetic field) can be found in Ref. [21]; here 
we just give its main steps. One starts with Eq. ((601. Wick's theorem gives 



(^(kK%K(sK™( P )) = 

(v i (k)v^(q))(v n ( S )v* m (p)) + 
(^(kK(s))(^'(q)u* m (p)) + 
(v\k)v* m (p))(v n (s)v*J( q )) . 



This has to be inserted in Eq. (60 ) together with the definition of the velocity power spectrum 



(54). Applying the projection operator gives the angular dependence 

abcd (k)[(6 a 



■q a q c )(5 bd - 
+ (S ad - q a q d ){5 hc - 
l + (k-(k^)) 2 + (k-q) 2 



(k - q) 6 (k - q) d ) 

(k^f(k^r)] = 

(k-q) 2 (k-(k^)) 2 . 



(112) 



Comparing then with Eq. (58), one arrives at the result: 



U v (k,t h t 2 ) = J d 3 pP v (p, t 1 ,t 2 )P u (|k-p|,t 1 ,i 2 )(l + 7 2 )(l + /3 2 ), (113) 



with 7 = k p, (3 = k k 
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